#### Pollution and Regime Support Analysis

rm(list = ls())

pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", "fracdiff",
         "fractal",
         "lme4", "ArfimaMLM","gam",
         "forecast",
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", 
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)

# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

# loading the data
load("ppp_cleaned")

# loading stock data
stock <- read.table(file = "TRD_index.txt", header = T)
# only focusing on 上证指数 
stock1 <- stock[stock$Indexcd == 1, ]
colnames(stock1)[2] <- "begin_date"
stock1$diff <- log(stock1$Clsindex) - log(lag(stock1$Clsindex))

merged <- merge(p2, stock1, by = c("begin_date"), all = T)
p2 <- merged[is.na(merged$aqi) == F, ]


for (i in 1:nrow(p2)) {
  if (is.na(p2$Clsindex[i]) == T) {
    p2$Clsindex[i] <- p2$Clsindex[i-1]
  }
}
p2$Clsindex


for (i in 1:nrow(p2)) {
  if (is.na(p2$diff[i]) == T) {
    p2$diff[i] <- p2$diff[i-1]
  }
}
p2$diff

p2$Clsindex <- (p2$Clsindex-mean(p2$Clsindex))/(sd(p2$Clsindex))
p2$diff <- (p2$diff-mean(p2$diff))/(sd(p2$diff))


##
# aqi and government satisfaction at the local level
l_satc <- lm(l_sat
             ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ  +
               + kids + married + income + gender + ccp + parade + national + age + Clsindex, data=p2)
m.vcov1 <-cluster.vcov(l_satc, p2$day)  
coeftest(l_satc, m.vcov1)
se11 <- coeftest(l_satc, m.vcov1)

NWvcov <- NeweyWest(l_satc, lag = NULL, prewhite = FALSE)
sen11 <- coeftest(l_satc, NWvcov) 

# aqi and c_sat
c_satc <- lm(c_sat
             ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ  +
               + kids + married + income + gender + ccp + parade + national + age + Clsindex, data=p2)
m.vcov25 <-cluster.vcov(c_satc, p2$day)  
coeftest(c_satc, m.vcov25) 
se22 <- coeftest(c_satc, m.vcov25)

NWvcov <- NeweyWest(c_satc, lag = NULL, prewhite = FALSE)
sen22 <- coeftest(c_satc, NWvcov) 

# aqi and the demand for more watchdog at the local level
l_watchc <- lm(l_watch
               ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ +
                 + kids + married + income + gender + ccp + parade + national + age + Clsindex, data=p2)
m.vcov13 <-cluster.vcov(l_watchc, p2$day)  
coeftest(l_watchc, m.vcov13) # coeftest
# robust p = .03313
se33 <- coeftest(l_watchc, m.vcov13)
NWvcov <- NeweyWest(l_watchc, lag = NULL, prewhite = FALSE)
sen33 <- coeftest(l_watchc, NWvcov) 

# aqi and the demand for more watchdog at the central level 
c_watchc <- lm(c_watch
               ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ +
                 + kids + married + income + gender + ccp + parade + national + age + Clsindex, data=p2)
m.vcov14 <-cluster.vcov(c_watchc, p2$day)  
coeftest(c_watchc, m.vcov14) # coeftest
# robust at .065
se44 <- coeftest(c_watchc, m.vcov14)
NWvcov <- NeweyWest(c_watchc, lag = NULL, prewhite = FALSE)
sen44 <- coeftest(c_watchc, NWvcov)

# aqi and anti-west sentiment
anti_westc <- lm(anti_west
                 ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ  +
                   + kids + married + income + gender + ccp + parade + national + age + Clsindex, data=p2)
m.vcov25 <-cluster.vcov(anti_westc, p2$day)  
coeftest(anti_westc, m.vcov25) 
se66 <- coeftest(anti_westc, m.vcov25) 

NWvcov <- NeweyWest(anti_westc, lag = NULL, prewhite = FALSE)
sen66 <- coeftest(anti_westc, NWvcov) 


###### ###### ####
### Making the CLS Index Level Table
###### ###### ####

column.labels <- c("Local Government", "Central Government", "Local Oversight", "Central Oversight", "Anti-West")

cov.labs.sci <-c("AQI", "AQI Lagged", "AQI Lagged 2","Hukou Status", "Regime Insider", "SOE Employee", "Education","Children","Married","Income", "Female", "CCP Member","Parade Period","Holiday","Age", "SCI Level")

cls.level <- list(l_satc,c_satc,l_watchc,c_watchc,anti_westc)
## TABLE A10
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/sci_index_level.tex")
stargazer(cls.level, font.size = "large", digits = 4, 
          column.labels = column.labels, 
          covariate.labels = cov.labs.sci,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, keep.stat = c("n"),
          omit ="Constant",
          star.cutoffs = c(.10, .05, .01), 
          se = list(se11[,2], se22[,2], se33[, 2], se44[, 2], se66[, 2]),
          header = F, notes = "robust standard errors clustered by day in parentheses", float = F)
sink()

################ #
###### Models w/ CLS Difference Var
################ #

# aqi and government satisfaction at the local level
l_satc <- lm(l_sat
             ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ
             + kids + married + income + gender + ccp + parade + national + age + diff, data=p2)
m.vcov1 <-cluster.vcov(l_satc, p2$day)  
coeftest(l_satc, m.vcov1) # coeftest (this will give you results using clustered standard errors)
se11 <- coeftest(l_satc, m.vcov1)
# significant
NWvcov <- NeweyWest(l_satc, lag = NULL, prewhite = FALSE)
sen11 <- coeftest(l_satc, NWvcov) 

# aqi and c_sat
c_satc <- lm(c_sat
             ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
             + kids + married + income + gender + ccp + parade + national + age + diff, data=p2)
m.vcov25 <-cluster.vcov(c_satc, p2$day)  
coeftest(c_satc, m.vcov25) 
se22 <- coeftest(c_satc, m.vcov25)
# robust at .066
NWvcov <- NeweyWest(c_satc, lag = NULL, prewhite = FALSE)
sen22 <- coeftest(c_satc, NWvcov) 

# aqi and the demand for more watchdog at the local level
l_watchc <- lm(l_watch
               ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ
               + kids + married + income + gender + ccp + parade + national + age + diff, data=p2)
m.vcov13 <-cluster.vcov(l_watchc, p2$day)  
coeftest(l_watchc, m.vcov13) # coeftest

se33 <- coeftest(l_watchc, m.vcov13)
NWvcov <- NeweyWest(l_watchc, lag = NULL, prewhite = FALSE)
sen33 <- coeftest(l_watchc, NWvcov) 

# aqi and the demand for more watchdog at the central level 
c_watchc <- lm(c_watch
               ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ
               + kids + married + income + gender + ccp + parade + national + age + diff, data=p2)
m.vcov14 <-cluster.vcov(c_watchc, p2$day)  
coeftest(c_watchc, m.vcov14) # coeftest

se44 <- coeftest(c_watchc, m.vcov14)
NWvcov <- NeweyWest(c_watchc, lag = NULL, prewhite = FALSE)
sen44 <- coeftest(c_watchc, NWvcov) 



# aqi and anti-west sentiment
anti_westc <- lm(anti_west
                 ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ
                 + kids + married + income + gender + ccp + parade + national + age + diff, data=p2)
m.vcov25 <-cluster.vcov(anti_westc, p2$day)  
coeftest(anti_westc, m.vcov25) 
se66 <- coeftest(anti_westc, m.vcov25) 
# robust at .0765
NWvcov <- NeweyWest(anti_westc, lag = NULL, prewhite = FALSE)
sen66 <- coeftest(anti_westc, NWvcov) 

###### ###### ####
### Making the CLS Change Table
###### ###### ####

cov.labs.diff <-c("AQI", "AQI Lagged", "AQI Lagged 2","Hukou Status", "Regime Insider", "SOE Employee", "Education","Children","Married","Income", "Female", "CCP Member","Parade Period","Holiday","Age", "SCI Change")

cls.change <- list(l_satc,c_satc,l_watchc,c_watchc,anti_westc)
## TABLE A11
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/sci_index_change.tex")
stargazer(cls.change, font.size = "large", digits = 4, 
          column.labels = column.labels, 
          covariate.labels = cov.labs.diff,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, keep.stat = c("n"),
          omit ="Constant",
          star.cutoffs = c(.10, .05, .01), 
          se = list(se11[,2], se22[,2], se33[, 2], se44[, 2], se66[, 2]),
          header = F, notes = "robust standard errors clustered by day in parentheses", float = F)
sink()
